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ABSTRACT 

We present the results of a series of one-dimensional N-body and hydrodynamical 
simulations which have been used for testing the different clustering properties of 
baryonic and dark matter in an expanding background. Initial Gaussian random den- 
sity perturbations with a power-law spectrum P{k) oc k"' are assumed. We analyse 
the distribution of density fluctuations and thermodynamical quantities for different 
spectral indices n and discuss the statistical properties of clustering in the correspond- 
ing simulations. At large scales the final distribution of the two components is very 
similar while at small scales the dark matter presents a lumpiness which is not found 
in the baryonic matter. The amplitude of density fluctuations in each component de- 
pends on the spectral index n and only for n — —1 the amplitude of baryonic density 
fluctuations is larger than that in the dark component. This result is also confirmed 
by the behaviour of the bias factor, defined as the ratio between the r.m.s of baryonic 
and dark matter fluctuations at different scales: while for n = 1, 3 it is always less 
than unity except at very large scales where it tends to one, for n = — 1 it is above 
1.4 at all scales. All simulations show also that there is not an exact correspondence 
between the positions of largest peaks in dark and baryonic components, as confirmed 
by a cross-correlation analysis. The final temperatures depend on the initial spectral 
index: the highest values are obtained for n = —1 and are in proximity of high density 
regions. 

Key v^rords: galaxies: clustering - galaxies: formation ~ large-scale structure of Uni- 
verse 



1 INTRODUCTION 

In recent years a considerable improyement in the study of 
large-scale structure has been obtained with the develop- 
ment of three-dimensional codes able to follow simultane- 
ously the evolution of dark and baryonic matter. It is im- 
portant to include the gas in N-body simulations for at least 
two reasons: first we must estimate the large amount of gas 
that is observed in clusters of galaxies through X-ray emis- 
sion; second, we must remove any ambiguity in the iden- 
tification of galaxies. In the numerical simulations, the hy- 
drodynamical eq uations have been solved either using mesh 
based methods (IChiang, & Vishniac 1989|; Cen 1992 



Ryu et al. 1993; |Bryan et al. 1994) or particle methods like 



the s i iioothed part i cle hydrodynamics ( Hernquist 
198S|; lEvrard 199(i[ 



Katz 



Katz fc Gunn 1991). The latter, being 



Lagrangian methods, have the advantage of covering a larger 
range of clustering scales, but they are not very reliable in 
low density regions for the calculation of thermodynamical 



quantities. Eulerian methods, on the contrary, give a better 
estimate of thermodynamical quantities, but are constrained 
by the limited grid resolution in their capability of follow- 
ing the collapse of high density regions. A large dynamical 
range can be attained only by increasing the number of grid- 
points, but this is obviously limited by the available com- 
puter memory. At present the largest mesh used in Eulerian 
calculations has 256^ grid-points. The analysis of the evolu- 
tion at different cosmological scales is obtained by changing 
the physical normalization of the box size, and sometimes 
the results obtained at a given scale are used for setting up 
initial and boundary conditions at another scale (e.g. Cen 
& Ostriker 1992). 

Since three-dimensional calculations pose severe lim- 
its on the number of grid-points used in Eulerian methods, 
there has been in recent years an attempt to make up for the 
lack of resolution with the use of more accurate algorithms. 
One is looking for numerical methods able to well describe 
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both low density regions and strong gradients in the fluid 
flow. Very interesting is, for example, the possibility of intro- 
ducing mesh refinement schemes in grid-based codes which, 
associated with shock capturing methods, would consider- 



ably i mprove the resolution in high density regions (Anni- 
nos, Ncrman & Clarke 1994). 



The limited resolution reachable in three-dimensional 
codes makes sometimes difficult to separate, in the discus- 
sion of results, those effects which are only numerical from 
those which contains real physical information. Analysis like 
the one of Kang et al. (1994) are very useful for discrimi- 
nating between different methods in three dimensions, but, 
because of the limitation in computer memory affecting all of 
the methods, a simultaneous one-dimensional analysis of the 
same phenomena appears quite convenient. Working only 
in a single dimension makes possible to span a large dy- 
namical range and to follow structure formation at different 
scales. Approximate analytic solutions and their behaviour 
in highly non-linear regimes can also be better tested in 
one-dimensional calculations (e.g. Williams et al. 1991). Re- 
cently, Quilis, Ibanez & Saez (1994) have studied the appli- 
cability of modern high-resolution shock-capturing methods 
to the study of cosmological structures in presence of pres- 
sure forces. However, they do not follow at the same time 
the dynamics of the dark matter and so they are not able 
to compare the clustering properties of the two different 
components. In our work, instead, we study the evolution 
of cosmological one-dimensional perturbations when both 
baryonic and dark matter are considered. 

A detailed study of the evolution of gas and collision- 
less matter in a single pancake was presented by Bond et 
al. (1984) and Shapiro & Struck-MarceU (1985). In both of 
these works one-dimensional calculations of the coupled evo- 
lution of baryonic and dark component were performed in- 
cluding the effects of radiative and Compton cooling terms 
and thermal conductivity of the gas. The main difference 
between these two analyses is in the choice of the numeri- 
cal algorithm used for solving the hydrodynamical equations 



whi ch w as Lagrangian in one case ( Shapiro fc Struck-Marcel] 
198£|)~a|nd Eulerian in the other one ( |Bond et al. 1984| ). A 
similar analysis has been recently done by Thoul & Weinberg 
(1994), for a spherical configuration. They also have used 
a Lagrangian approach which has the advantage of giving 
a higher resolution where needed without having, in one- 
dimension, the same problems of grid distortion present in 
more dimensions. All of the previous one-dimensional works 
following the coupled evolution of dark and baryonic mat- 
ter consider isolated perturbations in the expanding uni- 
verse and start their simulations when hydrodynamical ef- 
fects begin to become important in the evolution. The infor- 
mation concerning the formation epoch and mass of these 
isolated structures is derived from N-body calculations as- 
suming some model of large-scale structure formation. 

In this work we do not limit our attention to isolated 
perturbations but we want to explore the differences in the 
statistical properties of the one-dimensional clustering of 
baryonic and dark matter (hereafter BM and DM, respec- 
tively) components in an expanding background. We anal- 
yse the effect of pressure forces and adiabatic heating in the 
dynamics of the gas and we do not include, at this stage, 
cooling terms in our computations. Our aim is to under- 
stand first the evolution of one-dimensional structures only 



in the presence of "compressible" effects. In addition, the 
inclusion of cooling processes requires the specification of 
characteristic time and length scales in the problem and in 
one-dimensional cosmological calculations it is not obvious 
how one can choose these quantities in connection with the- 
ory and observations. Therefore, we have preferred to make 
simulations which are scale- free; the absence of features in 
the primordial power-spectrum together with the higher spa- 
tial resolution implied by the use of one-dimensional simu- 
lations would permit to better discriminate between numer- 
ical artifacts and physical effects in the results. At the end, 
physical information on the computed quantities can easily 
be extracted by fixing the normalization quantities, like the 
box size and the final time. 

The plan of the paper is as follows. In Section 2 we 
introduce the dynamical equations both for the hydrody- 
namical variables and for coUisionless matter. The numer- 
ical methods used for the solution of these equations are 
also presented together with some numerical tests. The re- 
sults obtained in numerical simulations of cosmological one- 
dimensional structures are shown in Section 3. Discussion 
and conclusions are drawn in Section 4. 



2 DYNAMICAL EQUATIONS 

In this section we present the dynamical equations that we 
use for evolving the BM and DM components in the case of 
one-dimensional perturbations. We assume a flat universe 
{Q, = 1) in which the BM component accounts for 10% of 
the total mass {^^^.j = 0.1) and the rest is in the form of 
DM (^Djvj = 0.9). In the present analysis we neglect the 
effects of radiative and Compton coolings and any possible 
external heating. The DM component is approximated as a 
pressureless fluid and the BM component as a perfect fluid. 



2.1 Hydrodynamical equations 

We introduce the following set of dimensionless variables: 



t 

X 

X = —, 
XO 

to 

V = V , 

XO 

Q = , 

Qo 

3 to 

e = a e 



Qoxl ' 



to 



(1) 
(2) 
(3) 
(4) 
(5) 
(6) 



r 2/3 



where x is the comoving coordinate, t is the time, a — t 
the cosmological expansion factor, </> is the peculiar gravita- 
tional potential, g and e are the matter and internal energy 
densities, respectively. We have three basic units of normal- 
ization: the final time io, the comoving cell size xo and the 
mean baryon density go at time to. For a cosmological ap- 
plication we can take to and go equal to the present age and 
density of the Universe, respectively. 
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The evolution of the coUisional component is described 
by the following set of hydrodynamical equations; 

dg 1 d 



dt a dx 

dgv 1 d |...2 -1 CL 1.90 
dt a dx ^ a ^ dx 



dg 1 d 



a 1 dv 
- -p— 
a a dx 



(7) 
(8) 
(9) 



where p is the comoving pressure. 

We use the equation of state of an ideal gas with adia- 
batic index 7 = 5/3, so that 



p = (7 - l)e 



(10) 



The physical temperature T is related to the previous 
normalized quantities by 



(11) 



where Kb is Boltzman's constant and the mass of the 
proton. 

2.2 Collisionless matter 

In the dynamical equations for the collisionless particles we 
use as time integration variable the scale factor o: 

dx 

da ' 

3 . 9 , , 

1 u — — ^ . (13) 

da 2a 4a dx 

The peculiar gravitational potential due to both the 
BM and DM component is computed by solving the Poisson 
equation: 



(12) 



dx^ 



K ,^ 

= —\Qbm + Qe 



1) 



(14) 



where K = A-KCgotf-^ and G is the Newton's constant. 
2.3 Numerical Methods 

The hydrodynamical equations have been int egrated using 
the F l ux Corrected Transport (FCT) method ( Boris & Book 
197S; Zalesak 1979), a hybrid shock capturing method. In 
this method two difference schemes are blended together: 
a second order Lax-Wendroff scheme is used in regions of 
smooth flow, while a first order Lax scheme is used near 
discontinuities. Then second order accuracy is ensured ev- 
erywhere except near flow jumps where the dissipation in- 
troduced by the low order scheme guarantees monotonicity 
in the behaviour of flow variables. The blending between the 
two schemes is controlled by a monotonicity constraint that 
leads to the sharpest possible discontinuity proflles. How- 
ever, if the grid is not sufficiently fine, very strong gradients 
can be represented as a sequence of discontinuous jumps. 
This effect usually does not hinder the convergence of the 
method a nd can be reduced or com pletely avoided refining 
the grid ( [Woodward fc Colella 198^ ). 



For the integration of the collisionless matter we have 
used a particle-mesh code (Hockney & Eastwood 1981). We 



code, while the number of particles is a multiple of the num- 
ber of mesh points in order to have a good resolution also 
in low density regions. The interpolation used for computing 
the mass density and the forces acting on each particle is ob- 
tained by a TSC scheme (e.g. Hockney & Eastwood 1981), 
which ensures a good accuracy in the estimates of previous 
quantities without leading to an excessive slowing down of 
the computation. For the time integration we use a second 
order leap-frog method so that the accuracy in the N-body is 
comparable to that attained in the hydrodynamical part of 
the code. The peculiar gravitational potential is computed 
using the standard fast Fourier transform technique. 

The accuracy of the code has been tested against known 
analytical solutions. For the hydrodynamical part we per- 
formed the shock tube test comparing the analytical and 
numerical solutions. We found a very good agreement in 
smooth regions, while the shock is numerically diffused over 
about five grid-points. The evolution of an initial sinusoidal 
perturbation has been instead used for testing the perfor- 
mance of the N-body code: in this case the numerical solu- 
tion has been compared wi th the results ob tained using the 
Zel'dovich approximation ( Zel'dovich 197[l( ) which provides 
the exact one- dimensional solution up to shell crossing. In 
Figure 1 we compare, just before shell crossing, the result- 
ing DM distributions for the Zel'dovich and N-body solu- 
tions for two different choices of the number of grid-points 
Ng {Ng = 2^° and Ng = 2^^). The number of particles is 
instead kept fixed and set equal to 2^®. We notice that al- 
ready with the smallest number of grid-points the N-body 
solution is in very good agreement with the Zel'dovich re- 
sult. The difference in the height of the peak in the two cases 
is only due to the different numerical resolution implied by 
the adopted grids. 

We used a sinusoidal perturbation also for testing the 
coupled evolution of dark and baryonic components and to 
see how far the code is able to follow the collapse of a bary- 
onic fiuctuation. In Figure 2 we show the density profile of 
both components after the shell crossing time, obtained us- 
ing a grid of 2^"' points. We can see that fiuctuations in the 
baryonic matter as larger as lO'^ can be easily resolved with 
this resolution. Some typical problems of the FCT method 
start to appear after the formation of strong shocks which 
are not well resolved using the previous grid, but they do 
not destroy the solution and a local spatial refinement would 
permit to continue the evolution. Since situations as this one 
are never found in the following cosmological simulations, we 
believe that a grid of 2^^ is sufficient for our purposes. 



use the same mesh as in the hydrodynamical part of the 



3 COSMOLOGICAL ONE-DIMENSIONAL 
STRUCTURES 

In this section we present the results of the evolution of one- 
dimensional perturbations in an expanding background. We 
analyse the distribution of BM and DM, paying particular 
attention to the clustering properties of the two components 
and to the effects due to their coupling. 



3.1 Numerical Simulations 

We consider a grid of comoving length L, whose physical 
size grows according to the expansion parameter a. We sub- 
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divide our grid into 2^^ nodes and adopt periodic boundary 
conditions. The number of particles used in the N-body part 
of the code is 2^^. 

The initial conditions for a given model are determined 
as follows. The DM component is perturbed according to a 
Gaussian random density field characterized by a power-law 
spectrum of the general form 

P{k) = Ak"e-x.p{-k'^Rj) , (15) 

where ^ is a normalization constant and n the spectral in- 
dex. The short wavelength cut-off at the scale Rf — 2 grid- 
points ensures that the results are not affected by the sam- 
pling of modes whose size is close to that of the resolution 
of the simulation. 

The normalization constant A is fixed requiring that 
the onc-dirncnsional variance at the final time to (corre- 
sponding to the expansion factor o = 1) for linearly evolved 
perturbations is equal to unity. In computing the variance 
we assume a filtering radius J?* = 10~^L so that 

1 

cr'^{R^) = - P{k)W'^{k,R^,) dk . (16) 
Jo 

The function W{k,R^) is a one-dimensional top-hat filter: 
Wik,R.) = ^^^. (17) 

Initially the DM particles are evolved according to the 
Zel'dovich (1970) algorithm which is known to give the ex- 
act solution up to the shell crossing. We start our numerical 
computation when one of the following requirements is satis- 
fied; i) the largest density fluctuation is equal to unity; ii) the 
first shell crossing has occurred: this typically corresponds 
to epochs with a < 10~^ for all models. 

The numerical evolution of BM starts when a = 10~^, 
i.e. just after the recombination epoch. Before this time we 
assume that the distribution of BM is constant, as it would 
be in comparison with DM distribution, because of the cou- 
pling between radiation and baryonic matter before hydro- 
gen recombination. The internal energy is fixed assuming 
that the cooling of the gas is due only to adiabatic expan- 
sion between the epoch of recombination and the initial time 
of the simulations. In the subsequent evolution the pertur- 
bation in the baryonic matter is induced by gravitational 
couphng with the DM component. 

As already pointed out, our simulations are completely 
scale-free: the final amplitude of the fluctuations depend 
only on the scale J?* relative to which the initial normal- 
ization is flxed. The final thermodynamical quantities are 
completely determined once the amplitude of initial pertur- 
bations and the initial value of e are given. For example, if 
we think of R* as corresponding of 8 {h is the Hubble 
constant in units of 100 km s~^ Mpc~^), where observation- 
ally the variance is found to be unity, the whole box would 
correspond to L = 800 Mpc (we use i?* — L/lOO) and 
then we would be looking the evolution at very large scales. 
In order to study the evolution on smaller scales but fix- 
ing the normalization of the initial spectrum at the same 
physical scale, we should use a larger value for in units 
of L and, consequently, we should start with larger density 
fluctuations. These simulations would produce larger flnal 
temperatures on these smaller scales as we expect. 



3.2 Results 

We consider three difi'erent models with primordial spectral 
index n = —1, 1, 3: in three dimensions, these would corre- 
spond to n = —3, —1, 1 (the so-called Harrison- Zel'dovich 
spectrum) respectively, covering in this way the range of 
values usually adopted for cosmological models. For each 
model we run three simulations with difi'erent realizations 
of the initial conditions in order to obtain more accurate es- 
timates of the relevant quantities. Here we consider various 
statistical tests, such as the distribution of density perturba- 
tions, peaks and thermodynamical quantities, their correla- 
tions and the density power-spectrum. The results referring 
to the DM component have been previously smoothed by 
a Gaussian filter with a radius of ten grid-points, roughly 
corresponding to twice the numerical spreading of the shock 
front in the hydrodynamical calculations. 

3.2.1 Distribution of Density Perturbations 

In Figures 3a and 3b we show the behaviour of and 
in a realization for each model at two diiferent epochs: 
a = 0.5 and the final one a = 1, respectively. The results 
refer to about a tenth of the whole grid. Different runs of 
the same model show a very similar qualitative behaviour. 
The density fiuctuations are computed with respect to the 
mean value of the corresponding component. As expected, 
increasing the spectral index n, we have less power on large 
scales and, consequently, peaks in density contrast are more 
frequent, but they have a smaller density contrast, while 
underdense regions are less extended. We notice that only 
for n = — 1 peaks in BM are higher than those in DM and 
their contribution to the local gravitational field starts to be 
comparable to that of DM. 

For the n = 3 case there is a correspondence between 
the distribution of DM and BM although there are always 
substructures in the coUisionless component which are not 
present in the other one (see also the following discussion on 
relative bias, power-spectrum and number of peaks). This 
behaviour is enhanced as we decrease the spectral index: 
BM appears more and more clumped while DM stays more 
spread with small substructures. This is particularly evident 
in the structures appearing in the model with n = —1, where 
coUisionless matter is spread over a region which is several 
times larger than that occupied by the peak in BM. It is 
interesting to notice also that there is not an exact corre- 
spondence between the positions of the largest BM and DM 
peaks. 

Comparing the figures at two different times we see the 
tendency to the merging of substructures which produces 
higher peaks and larger voids. This behaviour is again more 
evident in simulations with more power on large scales. For 
example, in the panels referring to the n = —1 simulation, 
it is possible to note that the two central peaks in the BM 
component at a = 0.5 quickly merge forming a unique high- 
density region at a = 1. This phcuomenon is less pronounced 
for DM, which at o = 1 continues to be divided into small 
substructures. 

A (luantitative estimate of the spatial correspondence 
between the BM and DM density distributions can be ob- 
tained studying the cross-correlation coefficient (see e.g. 
Coles, Melott & Shandarin 1993) 
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S{R) = 



(18) 



where 5„., and 5, 



are respectively the density contrasts 



in the BM and DM components, smoothed with a Gaussian 
window of radius R and ai = (<5|)^/^ are the corresponding 
r.m.s.. The mean is computed over all grid-points. Defini- 
tion (18) implies | 5 |< 1. The limit S = -1-1 corresponds to 
^BM ~ ^^DM' where C is a constant: in this case there is 
a perfect agreement in the positions of the structures in the 
two components. In Figure 4, the cross-correlation coefficient 
S is plotted as a function of the Gaussian filtering radius R 
for the three models when a — 0.5 and a = 1. The error bars, 
shown for clarity only at the final time, represent the scatter 
r.m.s. between the three different realizations of each model. 
At earlier times these errors are found to be always smaller. 
As expected, decreasing the filtering radius and/or allow- 
ing for a longer evolution, the differences between the two 
components increase and, consequently, the coefficient S de- 
creases. This is particularly true for the model with n = —1, 
where, for small R (up to L/400), the cross-correlation is 
even less than 0.4, denoting the small correspondence be- 
tween BM and DM distributions which was already clear 
in Figure 3. On the contrary, a value of S close to unity is 
obtained in the simulations with n = 3, even when a small 
filtering radius is adopted. 

The study of high-density regions in BM is particularly 
interesting because they are expected to be related to the po- 
sitions where galaxy formation occurs. We define as peaks all 
the grid-points which are local maxima of the density fluctu- 
ation field. In Figure 5 we show, at two different times, the 
number of peaks Npk having an height greater than a given 
threshold S, both for BM and DM components. The com- 
parison between the models with different spectral indices 
shows once again a different behaviour. Only for the model 
with n = — 1, the number of high peaks in BM is larger than 
that of DM at both considered times. The latter model is 
also the one which shows the largest time evolution in the 
number of peaks. In the other cases, Npk does not changes 
significantly and the peak number in the BM component is 
always less than that in DM one. 



3.2.2 Thermodynamical quantities 

In Figure 6 we show the behaviour of the baryonic quantities 
Q^^j , p and T for the case n = — lata = l (the results refer 
to the same realization shown in the upper panel of Figure 
3b but in this case the whole grid is displayed). This is the 
model in which the highest values in thermodynamical quan- 
tities are obtained, thus some effects, which are also present 
in the other models, are better visible. The baryonic density 
and the pressure are given in normalized units, while the 
temperature is shown in Kelvin degrees. In absence of heat- 
ing from external sources and cooling processes, wc notice 
a strong correspondence among peaks in density, pressure 
and temperature. However, in correspondence of the largest 
density fluctuations we observe a double peak on the top 
of the temperature profile which shows, as expected, that 
heating mainly occurs in the regions in which the velocity 
gradient is large. 

The temperature distribution in the different models 
is better illustrated by the cumulative volume filling factor 



and mass fraction F as a function of the temperature T. 
Figure 7 shows the behaviour of these quantities at a = 1 
as obtained by summing over the contribution of all three 
realizations for each model. We notice that the fraction of 
matter at high temperature (e.g. T > lO** K) increases as we 
decrease the spectral index: it is approximately 15% and 2% 
for n = — 1 and n = 1, respectively, while no grid-point has 
such a high temperature in the n = 3 model (in this case 
the highest temperature obtained is about 25 K). At the 
same time, most of the volume is at low temperatures: this 
is particularly evident in the case n = — 1, where, due to the 
dominance of empty regions, more than 90% of the volume 
has T < 10^^ K. On the contrary, the two cumulative func- 
tions are more similar in the model with n = 3, where struc- 
tures are smaller and more uniformly distributed. These re- 
sults suggest that, in absence of other dissipative or heating 
effects (viscosity, thermoconduction, starbursts, etc.), it is 
difficult to obtain from one-dimensional cosmological large- 
scale perturbations temperatures as high as those seen in 
X-ray observations. 

Figure 8 shows contour plots for the number of zones 
characterized by a given temperature and baryonic density. 
The additional information which we can extract from these 
graphs is the correlation between these two quantities in 
all the three models. This correlation is particularly strong 
in the case n = —1 and refers to the large empty regions 
characterizing this model. On the contrary, more scattered 
distributions are presented in the models with n = 1 and 
n = 3. 

3.2.3 Power-spectrum and bias 

In order to study the time evolution of clustering at differ- 
ent scales, we have computed the power-spectrum P{k). If 
the density fluctuations remain linear at large scales, P{k) is 
expected to grow as . In Figure 9 we show the behaviour 
of the power-spectrum both for DM and BM for the three 
models at different times. The units of the wavenumbers k 
are such that k = 1 corresponds to the fundamental wave- 
length of the computational grid. Considering the models 
with n — 3 and n = 1, we see that at small wavenumbers 
(i.e. at large scales) the growth of the perturbations in the 
DM component is in good agreement with linear theory. The 
situation is different for the model with n = — 1 where the 
growth of perturbations is non-linear even at these scales. 
At larger wavenumbers (i.e. at small scales) deviations from 
linear theory are found, as expected. In particular, we notice 
that there is a faster evolution for the models with n = — 1 
and n = 1: the non- linear growth produces the coupling of 
high- and low-fc modes. In the model with n = 3, where 
initially the large-scale power is small, this effect does not 
appear. Then, our simulations, even if starting with a power- 
law spectrum, seems to exclude the possibility of having self- 
similar evolution when large scale power is present. Only the 
model with n = 3 presents a self-similar evolution in the DM 
at least for log k < 2.4. 

The behaviour of the BM is very different. At the begin- 
ning, due to the assumption of uniform distribution on all 
scales, the power-spectrum is zero. At early times all mod- 
els show a very slow growth produced by the coupling with 
the DM fluctuations (the curves for the model n = — 1 are 
not visible because P{k) is too low). Only at late epochs 
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(a « 0.5), the fluctuations in the BM grow considerably and 
quickly reach those in the DM. At the final time a = 1 the 
power-spectra in the two components are different at small 
scales, where the growth of gas perturbations is delayed by 
the pressure: this is particularly true in the models with 
n = — 1 and n = 1, for which higher values of the pressure 
are reached. 

An interesting way to measure the relative growth of 
the perturbations is the bias factor 6, defined as 

b^R) = Z^M^ (19) 

where a^{R) is the variance of the density fluctuation field 
in the i-componciit, smoothed using a Gaussian window of 
radius R. Figure 10 shows the behaviour of h for the three 
models at a = 0.5 and a = 1. The scatter r.m.s. between 
the different realizations of each model is also plotted only 
at the final time. The models with n — 1 and n = 3 behave 
in a similar way: the DM appears more clustered than the 
BM at small scales, while at larger scales they tend to have 
the same clustering properties. As expected, the pressure in 
the BM prevents shell-crossing and a collapse comparable 
to that occurred in the DM. At larger scales where dissipa- 
tion is negligible, the two distributions are similar. In the 
case n = — 1 instead, the BM at the final time a = 1 has 
a larger variance than the DM component at all scales; this 
is in agreement with the density distribution shown in Fig- 
ure 3b. The motivation for the different behaviour of this 
model is the large peculiar velocities of the DM component. 
In fact, the absence of dissipative effects delays the forma- 
tion of more compact structures. In the BM, instead, higher 
temperatures and pressures are reached, but they are not 
sufficient to overcome the gravitational force and avoid the 
collapse. In addition the conversion of the kinetic into inter- 
nal energy allows the formation of more compact baryonic 
structures. 

From Figure 10 it is possible to extract also information 
about the time evolution of the bias factor. Again we can 
note a different behaviour between the models. In the case 
n = — 1, the evolution of b from a = 0.5 to a = 1 is rapid at 
all scales: b changes by a factor ~ 1.5. On the contrary, the 
increase of b in the same time interval is smaller for n = 1 
and n = 3 (approximately 15% and 5%, respectively) but 
the final value is always below unity. 



nent presents a lumpiness which is not found in the BM. A 
quite important feature shown by all simulations is the non- 
exact correspondence between largest DM and BM peaks. 

In our simulations peaks in the DM are always higher 
than those in BM except for n = —1, i.e. for the spectrum 
with more power on large scales. In this case the variance in 
baryon density distribution is laxger than that of DM at all 
scales. Moreover, in the highest baryon peaks the self-gravity 
of baryons is comparable to the gravitational field due to 
the coUisionless component and this can have an important 
effect on the evolution of the latter. In general, many peaks 
in the DM component are associated to a single peak in 
the BM. These are usually embedded in the corresponding 
baryonic peak which stays more spread in space because of 
the effect of pressure working against the collapse. For the 
highest BM peaks present in the case n = — 1 the situation is 
reversed: the pressure, although higher than in the previous 
cases, is not able to oppose sensibly the infall of BM and we 
end up with a distribution in the DM which is more spread 
than that in the BM. This behaviour is confirmed by the 
value of the bias which, for n = —1 is larger than one at all 
scales. 

Our analysis is completely scale-free and the final distri- 
bution of the gas temperature depends only on the assumed 
power spectrum. The highest values of the temperature are 
obtained for the case n = —1 and they correspond to the 
high density regions. As the spectral index increases, large- 
scale motion is reduced and this decreases the height of the 
baryonic peaks and therefore also the final heating of the gas 
is smaller. Dissipative processes present in the coUisional 
component allow the transformation of kinetic into inter- 
nal energy and this makes possible the formation of narrow 
baryonic peaks even in the absence of cooling processes not 
included in the present work. In our simulations pressure 
forces have never been able to equal gravitational forces and 
so the motion of infall is expected to contirme even in pres- 
ence of high temperatures. The inclusion of cooling terms 
would amplify this effect and possibly cause a fragmentation 
of baryonic peaks. Work is presently in progress for testing 
the effect of cooling and heating processes in the evolution of 
one-dimensional cosmological perturbations. This, however, 
will require the specification of the physical scales involved 
in the problem and therefore the scale-free character of the 
present study will then be lost. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have studied the evolution of one- 
dimensional perturbations in a medium composed by BM 
and DM. Initial Gaussian perturbations with a power-law 
spectrum are assumed in the dominant DM component and 
then the perturbations in the BM component, initially uni- 
formly distributed, are induced through gravitational cou- 
pling. We observe that BM does not fall immediately into 
the potential wells created by the distribution of the coUi- 
sionless component, but when this happens, the amplitude of 
density perturbations in the gas becomes quickly compara- 
ble to that of the DM. At large scales the flnal distribution 
of the two components, as shown by the cross-correlation 
coefficient and by the bias factor, is very similar. The main 
differences are present at small scales where the DM compo- 
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at the time a = 1 for the different models: n — ~1 (left), 
n = 1 (centre) and n = 3 (right). The contour levels are 
defined as follows: lO'^'''''', where / is a positive integer; the 
outermost contour corresponds to / = 4 and contours inside 
it have gradually increasing I. 

Figure 9. The time evolution of the power-spectrum P{k) 
for baryonic (left column) and dark matter (right column) 
component for the different models: n = — 1 (top), n = 1 
(centre) and n — 3 (down). The curves refer to different 
epochs: a = 0.09, 0.17, 0.25, 0.33, 0.5 (dotted lines from down 
to top) and a = 1 (solid line). 

Figure 10. The relative bias b = (^bm I'^dm ^ ^ function of 
the filtering scale R aX a = 0.5 (dotted line) and a = 1 (solid 
line) for the different models: n = ~1 (left), n — 1 (centre) 
and n = 3 (right). Error bars (displayed for clarity only 
at a = 1) represent the r.m.s. scatter between the different 
simulations. 



FIGURE CAPTIONS 

Figure 1. Sinusoidal test for the single dark matter compo- 
nent: distribution of density fluctuations just before shell 
crossing. Comparison between the results obtained from 
N-body (solid line) and Zel'dovich approximation (open 
squares), when two different numbers of grid-points Ng are 
adopted: Ng = 2^" (left) and Ng = 2^^ (right). 
Figure 2. Sinusoidal test for the coupled evolution of bary- 
onic and dark matter components: distribution of density 
fluctuations after shell crossing. The dotted and solid lines 
refer to the baryonic and dark matter components, respec- 
tively. 

Figure 3a. The density fluctuation distribution at a = 0.5 
for the three different models: n = — 1 (top panel), n = 1 
(central panel) and n = 3 (bottom panel). The dotted and 
solid lines refer to the baryonic and dark matter compo- 
nents, respectively. Only about a tenth of the whole grid is 
displayed. 

Figure 3b. As Figure 3a, but at the final time a — 1. 
Figure 4. The cross-correlation coefficient S as a function 
of the Gaussian filtering radius R (in units of the whole 
grid L) at the time a = 0.5 (dotted line) and a = 1 (solid 
line) for the different models: n = —1 (left), n — 1 (centre) 
and n = 3 (right). Error bars (displayed for clarity only 
at a = 1) represent the r.m.s. scatter between the different 
simulations. 

Figure 5. The number of peaks Npk in dark (solid line) 
and baryonic (dotted line) components as a function of their 
height S when a = 0.5 (left column) and a = 1 (right col- 
umn) for the different models: n = — 1 (top), n — 1 (centre) 
and n = 3 (down). 

Figure 6. The distribution of baryonic density fluctuation 
^BM' pressure p and temperature T for the case n = —1 at 
a = 1. The whole grid is displayed. Pressure is in arbitrary 
units while temperature is in Kelvin degrees. 
Figure 7. The cumulative volume filling factor (solid line) 
and mass fraction (dotted line) as a function of the tem- 
perature T at the final time a — 1 for the different models: 
n = — 1 (left), n — 1 (centre) and n — 3 (right). 
Figure 8. Contour plots for the number of grid-points char- 
acterized by a given temperature T and baryonic density q 
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